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3.1 



O Abstract 

Using direct numerical simulations (DNS) and mean-field simulations (MFS), the effects of non-uniformity of the magnetic field 
jEf^on the suppression of the turbulent pressure is investigated. This suppression of turbulent pressure can lead to an instability which, 
-in turn, makes the mean magnetic field even more non-uniform. This large-scale instability is caused by a resulting negative 
^< "contribution of small-scale turbulence to the effective (mean-field) magnetic pressure. We show that enhanced mean current density 
increases the suppression of the turbulent pressure. The instability leads to magnetic flux concentrations in which the normalized 
effective mean-field pressure is reduced to a certain value at all depths within a structure. 

Keywords: Pacs: 91.25.Cw, 92.60.hk, 94.05.Lk, 96.50.Tf, 96.60.qd 
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1. Introduction 

The Sun's magnetic field is generally believed to be due to a 
turbulent dynamo operating in the convection zone in its outer 
30% by radius [Brandenburg and Subramanian, 2005, Krause 
and Radler, 1980, Moffatt, 1978, Parker, 1979, Zeldovichef al., 
1983]. Recent simulations performed by a number of groups 
have shown that the magnetic field is produced in the bulk of 
the convection zone. According to the flux tube scenario, most 
of the toroidal magnetic field resides at the bottom of the con- 
vection zone, or possibly just beneath it [Gilman and Dikpati, 
2000, Parfrey and Menou, 2007]. Another possibility is that 
most of the toroidal field resides in the bulk of the convec- 
tion zone, but that its spatio-temporal properties are strongly 
affected by the near-surface dynamics Kapyla et al. [2012a], or 
the near-surface shear layer [Brandenburg, 2005]. In any case, 
the question then emerges how one can explain the formation of 
active regions out of which sunspots develop during the lifetime 
of an active region. 

In the past, this question was conveniently bypassed by re- 
ferring to the possible presence of a strong toroidal flux belt at 
the bottom of the convection zone, where they would be in a 
stable state, except that every now and then they would become 
unstable, for example to the clamshell or tipping instabilities 
[Cally et al., 2003]. However, if the magnetic field is contin- 
uously being destroyed and regenerated by the turbulence in 
the convection zone proper, the mechanism for producing ac- 
tive regions and eventually sunspots must be one that is able 
to operate within a turbulent environment. One such mech- 
anism may be the negative effective magnetic pressure insta- 
bility (NEMPI) which is based on the suppression of turbulent 



pressure by a weak mean magnetic field, leading therefore to 
a negative effective (or mean-field) magnetic pressure and, un- 
der suitable conditions, to an instability [Kleeorin et al, 1989, 
1990, 1996, Kleeorin and Rogachevskii, 1994, Rogachevskii 
& Kleeorin, 2007]. This has been the subject of intensive re- 
search in recent years [Brandenburg et al, 2010, 2012, Kapyla 
et al., 2012b, Kernel et al., 2012a,b,c, Losada, 2012], follow- 
ing the first detection of such an instability in direct numeri- 
cal simulations [DNS; see Brandenburg et al, 2011]. Another 
mechanism that has been discussed in connection with the pro- 
duction of magnetic flux concentrations is related to the sup- 
pression of the turbulent convective heat flux [Kitchatinov & 
Mazur, 2000]. Meanwhile, simulations of realistic solar con- 
vection have demonstrated that large-scale magnetic flux inho- 
mogeneities can develop when horizontal magnetic flux is in- 
jected at the bottom of the simulation domain [Stein and Nord- 
lund, 2012], but it remains to be seen whether this is connected 
with any of the two aforementioned mechanisms. 

The purpose of the present paper is to investigate the pos- 
sibility that higher-order contributions (involving higher spa- 
tial derivatives of the mean magnetic field) might play a role 
in NEMPI. We do this by using DNS to measure the result- 
ing turbulent transport coefficients in cases where a measur- 
able mean current density develops in the DNS. Note that even 
for an initially uniform mean magnetic field, a mean current 
density develops as a consequence of NEMPI itself, which re- 
distributes an initially uniform magnetic field into a structured 
one. To investigate this process further, we use appropriately 
tailored mean-field simulations (MFS) that show how the spa- 
tial variations of the negative effective magnetic pressure vary 
in space as the instability runs further into saturation. We be- 
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gin by discussing first the basic equations and turn then to the 
results. Throughout this work we use an isothermal equation 
of state which yields the simplest possible system to investigate 
this process. 

2. Basic equations 

In this paper we use both DNS and MFS to study the effects 
of nonuniformity of the magnetic field of the development of 
NEMPI. In the DNS, the solutions turn out to have a large-scale 
two-dimensional pattern which can best be isolated using av- 
eraging over the y direction. Furthermore, by imposing a uni- 
form magnetic field, we can determine some of the turbulent 
transport coefficients that characterize the dependence of the 
Reynolds and Maxwell stress on the mean field. This is gener- 
ally done by determining the total mean stress 



n, 



■ pUiUj+ ±6ijB 2 -BiBj. 



(1) 



Here, the vacuum permeability is set to unity and overbars in- 
dicate y averages. This total mean stress has contributions from 
the fluctuations, 



n ;/ = puiUj + \6ijb 2 - btbj, 



(2) 



where u — U - U and b = B - B are the departures from 
the averaged fields. Here U and B are the mean velocity and 
magnetic fields, and 75 is the mean fluid pressure. This, together 
with the contribution from the mean field, namely 



n" = p Ui Uj + Stj (p+\B 2 )- BtBj - 2vp Sy, 



(3) 



yields the total mean stress tensor, i.e., Ily = n, ; +n, ; . The term 

Hjj depends only on the mean field and is therefore directly 
— f 

obtained in MFS, while ITy is caused by the fluctuating velocity 
and magnetic fields and requires a parameterization. It has a 

contribution independent of the mean field, n ; '- , and one that 

depends on it, AHy(B). Much of the recent work in this field 
focussed on the parameterization 



\q v (B) Sij B + q s (B) B t Bj + q g (B) B 2 g igj , (4) 



where g; = gjg is the unit vector in the direction of gravity. 
This difference in the mean stress, All,- -(B), is caused solely by 
the presence of the mean magnetic field B, where q p (J3) is found 
to be a positive function of B = \B\/B eq only, and, for weak 
magnetic fields, q p (B) is well in excess of unity for ReM » 
1, thus overcoming the magnetic pressure from the mean field 
itself. However the functions q s (B) and q g (B) were found to 
be small for isothermal turbulence. The net result for the sum 

By + Ally IS 



-f,0 



n, 7 * Uij + 6 uPe ff(B/B eq ) - BiBj, 



(5) 



where p e g = \ \l — q p (B)] Br is the mean effective magnetic 
pressure that is negative for/? < B cr it, where B cr i t ~ 0.5 (depend- 
ing on other details of the system). This results in a large-scale 



instability (NEMPI) and the formation of large-scale inhomo- 
geneous magnetic structures. 

In the nonlinear stage of NEMPI, the mean magnetic field 
becomes strongly nonuniform. This implies that the Maxwell- 
Reynolds stress tensor AYl^ may depend also on spatial deriva- 
tives of the mean magnetic field, i.e., 

f 

Ally = -\Stj q v B 2 + q s B,B } + q g B 2 gi gj + CiB Um B j<m 
+C2B m .iB m j + C3 (Bi m B„ h j + Bj m B„ U j), (6) 

where B/j = V;B,-. We decompose By into symmetric and 
antisymmetric parts: 



B U = (3B\ 



±£ 7 

2°/ Jin J in •> 



(7) 



where (3B)y = j(Bij + Bjj). Substituting equation (7) into 
equation (6) we obtain: 



Ally 



-\Sij q p B 2 + q s B^Bj + q g B 2 g igj 
-q s (j 2 d u - JiJA - qp (dB) im (dB) mj 
-qi (simi(dB) mj + s jm i(dB) mi jli. 



(8) 



Let us consider a mean magnetic field of the form B = 
(0,B v (x,z),0), so (dB)„ = JJ2, (dB) yz = -JJ2 and 



(dB) im (dB) 



"1 .1 



_0_ 

I -JJ Z 



-J X J ^ 

f 
-2 







(9) 



[s im i(dB) mj + e Jm i(dB) m j J, 

t 
-1 



-J X J Z 

Equations (8) and (9) yield: 



-Jxh \ 


t 



AH„ 

ah; 

\TT' 

Alt 



= -qi J z - jq p B , 

= qil 2 + (q s - j^p) B~ , 

= -cj]ll + {q g - \q p ) B' , 

= qjJ x J z , 



(10) 



(11) 
(12) 

(13) 
(14) 



2 —2 2 —2 

where q 3 (B ,J ) = q,+qi/2+q F /4, q\{B ,J) = qi/2-q 3 -q F /4, 

—2 —2 —2 —2 —2 —2 

q p = q p (B ,J),q s = q s (B , J ) and q g = q g (B , J ). Unfortu- 
nately, we have only 4 equations, but 5 unknowns, so we cannot 
obtain all the required transport coefficients independently. In 
the following, we can only draw some limited conclusions that 

will allow us to motivate a numerical assessment of the nonlin- 
2 2 

ear (B ,J ) dependence of q p . 



3. Results 

3.1. DNS model 

Following the earlier work of Brandenburg et al. [201 1] and 
Kernel et al. [2012b,c], we solve the isothermal equations for 



2 



the velocity, U, the magnetic vector potential, A, and the den- 
sity, p 

PT? = - c 'Vp + / XB +p(f + g) + V- (2v P S), (15) 



D/ 
8A 
dt 
dp 

dt 



UxB + rpl 2 A, 
-V-pU, 



(16) 
(17) 



where v is the kinematic viscosity, rj is the magnetic diffusiv- 
ity due to Spitzer conductivity of the plasma, B = Bo + V X A 
is the magnetic field, Bo = (0, Bo, 0) is the imposed uniform 
field, J = V X B/fj.o is the current density, po is the vacuum 
permeability, S,j = h(Uij + Up) - T<5yV • U is the traceless 
rate-of-strain tensor. The forcing function, /, consists of ran- 
dom, white-in-time, plane, non-polarized waves with a certain 
average wavenumber, kf. The turbulent rms velocity is approx- 
imately independent of z with u Tms = (h 2 ) 1 ^ 2 w 0.1 c s , where 
c s = const is the isothermal sound speed. The gravitational 
acceleration, g = (0, 0, —g) is chosen such that k\H p — 1, so 
the density contrast between bottom and top is exp(27r) « 535. 
Here, H p = c 2 lg is the density scale height and k\ = 2n/L is 
the smallest wavenumber that fits into the cubic domain of size 
L? . We consider a domain of size L x x L y x L z in Cartesian co- 
ordinates (x, y, z), with periodic boundary conditions in the x- 
and y-directions and stress-free, perfectly conducting bound- 
aries at the top and bottom (z - +L z /2). In the following 
we refer to kf/k\ as the scale separation ratio, for which we 
choose the value 30 in all cases. For the fluid Reynolds num- 
ber we take Re = u^/vkf = 18, and for the magnetic Prandtl 
number PrM = v/77 = 0.5. The magnetic Reynolds number is 
ReM = PrMRe. In our units, yuo = 1 and c s = 1 . The simulations 
are performed with the Pencil Code 1 which uses sixth-order ex- 
plicit finite differences in space and a third-order accurate time 
stepping method. We use a numerical resolution of 256 3 mesh 
points. In the MFS we also use 128 meshpoints, but because the 
MFS are two-dimensional, our resolution is 128 2 mesh points. 
In the model presented below, the z extent is however slightly 
bigger: z/H p is 8 instead of 2n. 

3.2. DNS results 

In all cases, we consider a weak imposed magnetic field in 
the y direction. In Fig. 1 we show y-averaged visualizations of 
the normal component of the magnetic field, B Y , together with 
the two components of the mean current density, J x = -8B y /dz, 
and J. = dBy/dx, and the normalized effective magnetic pres- 
sure, P eff = AfiL/Beq- 

As a result of NEMPI, the field in the xz plane gets concen- 
trated in one position and diluted in another; see the first panel 
of Fig. 1. This leads to an equilibrium in which the resulting 
reduction of the turbulent pressure is offset by a corresponding 
increase in the gas pressure and therefore a corresponding in- 
crease in the density. The nonuniformity of the magnetic field 
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Figure 1: Representations of By, P e ff, and the two components of / in the xz 
plane from a DNS with ReM = 18. 
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Figure 2: The 3 diagonal components of the tensor An,.,, (solid, black), ATI 

(dashed, blue), and AIL. (dotted, red), normalized by B^„ with a zero current 
density. 

implies a non-vanishing current density which is best seen in J z 
(lower right panel of Fig. 1), but this is mainly because of en- 
hanced fluctuations resulting from variations in the z direction. 

In Fig. 2 we show three diagonal components of the tensor 
— f — f — f , 
AII fl ., ATL, and All.., normalized by B^ q in turbulence with a 

zero mean current density. Fig. 2 demonstrates that the tensor 
— f 

ATLy in the direction of the mean magnetic field is different 
— f — f 

from the tensors AII rt and ALT in the directions perpendicular 
to the mean magnetic field. 

In the following we want to study the possible effects of cur- 
rent density on the resulting mean-field (or effective) magnetic 
— f 

pressure. We find that All,., vanishes, which implies that qj = 0; 
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Figure 3: Normalized diagonal components An f rr and AIlL as a function of 
8, for vanishing mean current density (black, solid line), for points with low 
(dashed, blue, 0.25 > J H^/B^ > 0.1) and higher current densities (dotted, red, 
0.25 i J 2 /H~B 2 Q ). 
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Figure 4: Like Fig. 3, but compensated by 1/(1 + J~ /k*B] q ) favkjH p = 4. 
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see Eq. (14). On the other hand, since All » ALT , the coef- 
ficient q g vanishes or is small. In Fig. 3 we show two diago- 

— f — f , 

nal components AII tv and ATE,, normalized by B~ q for differ- 
ent mean current densities. Inspection of Fig. 3 shows that the 

mean current density increases the negative minimum of the ef- 

— f 

fective magnetic pressure characterized by AII A . r This implies 
that an enhanced current density increases the effect of negative 
effective magnetic pressure, i.e., they intensify the formation 
of magnetic structures. Therefore, Eq. (11) allows us to deter- 
mine q p = -2An tt /ZT and the mean effective magnetic pres- 
sure p e ff = j \l - g p (/?)] B = AYl XJC + jB . In agreement with 
earlier studies [Brandenburg et al., 2012], we find a clear neg- 
ative minimum in P^iB) atB^ 0.25. However, as the current 
density increases, the minimum of !P e ff(jS) deepens, suggesting 
that NEMPI might turn out to be stronger than originally antici- 
pated based on the dependence !P e ff (B) that does not distinguish 
between strong and weak current densities. A reasonable fit to 
such a behavior would be of the form 



n„ = -£(i+/7*3^) 9 p(0)fl 



(18) 



where kj is a free parameter. In Fig. 4 we show that we can get a 
good fit to the data for k jH p = 4. Note further that equation (12) 
has two unknowns q\ and q s which cannot be determined for 



Figure 5: Dependence ofP^g(J3) on the value of g (in units of kic 1 ). Note that 
the depth of the minimum decreases with increasing gravity. 



such a simple configuration of the mean magnetic field, B = 
(0,B y (x,z),0). 

3.3. Gravity effects in DNS 

We mention in passing the effect of changing gravity. It is 
clear that increasing gravity enhanced the anisotropy of the tur- 
bulence, which seems to have a reducing effect on the negative 
effective magnetic pressure; see Fig. 5. The reason for this is 
at the moment not well understood. We emphasize that this ef- 
fect is connected with q p and not with q g that was introduced in 
equation (4). 

3.4. MFS 

Next, let us investigate the spatial variations of P e g in a cor- 
responding MFS. We use the parameters yS+ and B p that are ap- 
propriate in the regime investigated in the DNS above, namely 
ihsluansHp = 10~ 2 , corresponding to kfH p a 30, Bo/B eq o = 0.4, 
B+ = 0.32 and B p = 0.05. The result is shown in Fig. 6, where 
we plot in the upper panel the xz dependence of !P e ff. Since 
the domain is periodic in the x direction, we were able to shift 
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Figure 6: Representation of !P c ff in the plane using a MFS with ;/to / H lms H p = 
10~ 2 , corresponding to k s H p a 30, B /B cq0 = 0.4,/?* = 0.32 andyS p = 0.05. 

the position of the minimum such that it lies approximately at 
x — 0. The white vertical line near x — and the horizontal 
white line near z/H p = -4.3 indicate positions along which we 
plot in the next two panels P s g at three different times. 

Initially, the minimum of P e s(B) occurs at the height z/H p « 
-2.5, but at later times the minimum broadens and we have 
"Pes ~ -0.035 in the range -5.5 < z/H p < 2.5. In the last panel 
of Fig. 6 we show that the horizontal extent of the structure 
becomes narrower and more concentrated as time goes on. 



4. Conclusions 

The present results have shown that NEMPI tends to develop 
sharp structures in the course of its nonlinear evolution. This 
becomes particularly clear from the MFS presented in §3.4. The 
results of §3.2 suggest that this might have consequences of an 
intensification of NEMPI with increasing \J\, as was demon- 
strated using DNS. At present it is not clear what is the appro- 
priate parameterization of this effect. One possibility is that the 
J dependence enters in the same way as the B dependence, i.e., 
Yl xx = + J /k^Bl q )q p (J3) E~ , where we treat kj as a free 

parameter, although this might be a naive expectation given the 
small number of data points and experiments performed. 

The present results are just a first attempt in going beyond 
the simple representation of the turbulent stress in terms if the 
mean field along. Other important terms include combinations 
with gravity as well as anisotropies of the form J t Jj. Further- 
more, if there is helicity, one could construct contributions to 
the stress tensor using products of the pseudo-tensors JjBj and 
J jBj with the kinetic or magnetic helicity. Such a construc- 
tion obeys the fact that the Reynolds and Maxwell tensors are 
proper tensors. Such pseudo-tensors might play a role in the so- 
lar dynamo where the a effect is believed to play an important 
role. However, nothing is known about the importance or the 
sign of such effects. It would thus be desirable to have an accu- 
rate method that allows one to determine the relevant turbulent 
transport coefficients. 
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